library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(coda)
Rcpp::sourceCpp("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.cpp")
source("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.R")
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
## 
## Attaching package: 'reshape'
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
## The following object is masked from 'package:dplyr':
## 
##     rename
## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths
## Registered S3 method overwritten by 'pryr':
##   method      from
##   print.bytes Rcpp
## 
## Attaching package: 'pryr'
## The following objects are masked from 'package:purrr':
## 
##     compose, partial
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2021 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor

Data

6 distribution osservate, 4 gruppi

u  <- 30
x0 <- replicate(1, c(rnorm(u,3),rnorm( 3*u,10,.5)))
x1 <- replicate(1, c(rnorm(2*u,3),rnorm( 2*u,10,.5)))
x2 <- replicate(3, c(rnorm(2*u,-5),rnorm(2*u,10)))
x3 <- c(rnorm(2*u,3),rnorm( 2*u,10,.5),rnorm( 2*u,-5,.5),rnorm( 2*u,0,.5))

yd <- (c(x0,x1,x2,x3))
yg <- rep(1:6,c(rep(4*u,5),8*u))
plot(yd,col=yg)

nb <- 10000

Run the model

fixed concentration

m1 <- CAM(y_obser = yd,
            y_group = yg,
            K0 = 20,
            L0 = 20,
            prior = list(
              # hyperparameters NIG
              m0= mean(yd), 
              k0=1/var(yd), 
              a0=3, b0=1,
              # hyperparameters alpha and beta
              a_alpha=3, b_alpha =3,
              a_beta =3, b_beta = 3),
            nsim = nb,
            burn_in = nb,
            thinning = 1,
            verbose = F,
            fixedAB = T,
            kappa=0.5,
            cheap=F,
            seed=12345,
           sigma.prop.beta = 1)

MH marginal beta

m2 <- CAM(y_obser = yd,
            y_group = yg,
            K0 = 20,
            L0 = 20,
            prior = list(
              # hyperparameters NIG
              m0= mean(yd), 
              k0=1/var(yd), 
              a0=3, b0=1,
              # hyperparameters alpha and beta
              a_alpha=3, b_alpha =3,
              a_beta =3, b_beta = 3),
            nsim = nb,
            burn_in = nb,
            thinning = 1,
            verbose = F,
            fixedAB = F,
            kappa=0.5,
            cheap=F,
            seed=123456,sigma.prop.beta = 1,
            conditional.beta = F
  )

Gibbs conditional beta

m3 <- CAM(y_obser = yd,
            y_group = yg,
            K0 = 20,
            L0 = 20,
            prior = list(
              # hyperparameters NIG
              m0= mean(yd), 
              k0=1/var(yd), 
              a0=3, b0=1,
              # hyperparameters alpha and beta
              a_alpha=3, b_alpha =3,
              a_beta =3, b_beta = 3),
            nsim = nb,
            burn_in = nb,
            thinning = 1,
            verbose = F,
            fixedAB = F,
            kappa=0.5,
            cheap=F,
            seed=123456,sigma.prop.beta = 1,
            conditional.beta = T
  )

Alpha e Beta

plot(as.mcmc(m1$A_DP),main="M1") 

plot(as.mcmc(m2$A_DP),main="M2") 

plot(as.mcmc(m3$A_DP),main="M3") 

par(mfrow=c(1,2))
acf(m2$A_DP)
pacf(m3$A_DP)

par(mfrow=c(3,1))
plot(as.mcmc(m1$B_DP),main="M1") 

plot(as.mcmc(m2$B_DP),main="M2") 

plot(as.mcmc(m3$B_DP),main="M3") 

par(mfrow=c(1,2))
acf( m2$B_DP)
pacf(m3$B_DP)

Clustering

par(mfrow=c(1,3))
matplot(m1$NN.cz,type="l")
matplot(m2$NN.cz,type="l")
matplot(m3$NN.cz,type="l")

par(mfrow=c(1,1))
plot(density(m1$NN.cz[,1],type="l") ,col=1,lwd=2)
## Warning: In density.default(m1$NN.cz[, 1], type = "l") :
##  extra argument 'type' will be disregarded
lines(density(m2$NN.cz[,1],type="l"),col=2,lwd=2)
## Warning: In density.default(m2$NN.cz[, 1], type = "l") :
##  extra argument 'type' will be disregarded
lines(density(m3$NN.cz[,1],type="l"),col=3,lwd=2)
## Warning: In density.default(m3$NN.cz[, 1], type = "l") :
##  extra argument 'type' will be disregarded

plot(density(m1$NN.cz[,2],type="l") ,col=1,lwd=2) 
## Warning: In density.default(m1$NN.cz[, 2], type = "l") :
##  extra argument 'type' will be disregarded
lines(density(m2$NN.cz[,2],type="l"),col=2,lwd=2)
## Warning: In density.default(m2$NN.cz[, 2], type = "l") :
##  extra argument 'type' will be disregarded
lines(density(m3$NN.cz[,2],type="l"),col=3,lwd=2)
## Warning: In density.default(m3$NN.cz[, 2], type = "l") :
##  extra argument 'type' will be disregarded

Estimated partition

r1=PSM(m1$Z_j)
image(r1)

a1=mcclust.ext::minVI(r1,method = "greedy")
plot(yd,col=rep(a1$cl,table(yg)))

r2=PSM(m2$Z_j)
image(r2)

a2=mcclust.ext::minVI(r2,method = "greedy")
plot(yd,col=rep(a2$cl,table(yg)))

r3=PSM(m3$Z_j)
image(r3)

a3=mcclust.ext::minVI(r3,method = "greedy")
plot(yd,col=rep(a3$cl,table(yg)))

Densities

library(reshape2)
## 
## Attaching package: 'reshape2'
## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast
## The following object is masked from 'package:tidyr':
## 
##     smiths
xx <- seq(-10,15,by = .1)
LLL <- list()

LLL1 <- posterior.densities(m1,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D    <- melt(LLL1)
p1 <- ggplot(D)+
  geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
  theme_bw()+scale_color_viridis_d()
LLL2 <- posterior.densities(m2,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D <- melt(LLL2)
p2 <- ggplot(D)+
  geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
  theme_bw()+scale_color_viridis_d()
LLL3 <- posterior.densities(m3,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D <- melt(LLL3)
p3 <- ggplot(D)+
  geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
  theme_bw()+scale_color_viridis_d()
library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
p1/p2/p3

par(mfrow=c(3,1))
r1 <- apply(m1$post.dens,c(1,2),mean)
matplot(r1,type="l",main="M1")
r2 <- apply(m2$post.dens,c(1,2),mean)
matplot(r2,type="l",main="M2")
r3 <- apply(m3$post.dens,c(1,2),mean)
matplot(r3,type="l",main="M3")